Indicators to assess temporal variability in marine connectivity processes: A semi-theoretical approach

Oceanographic connectivity in an effective network of protected areas is crucial for restoring and stabilising marine populations. However, temporal variability in connectivity is rarely considered as a criterion in designing and evaluating marine conservation planning. In this study, indicators were defined to characterise the temporal variability in occurrence, flux, and frequency of connectivity in a northwestern Mediterranean Sea area. Indicators were tested on semi-theoretically-estimated connections provided by the runs of a passive particle transport model in a climatological year and in three years between 2006–2020, showing large deviation from the climatological year. The indicators allowed comparing the temporal variability in connectivity of four zones, highlighted differences in connectivity due to their locations and the mesoscale hydrodynamics, and identified areas that require further investigation. The three indicators also showed that the temporal variability in connectivity was influenced by the duration and depth of particle transport, although no consistent pattern was observed in the indicator variations of the compared zones. Provided that specific objectives will be given when parameterising transport models (i.e., selection of focus species and time period), indicators of temporal variability in connectivity have potential to support spatial conservation planning, prioritise the protection of marine resources, and measure the effectiveness of Marine Protected Areas, in line with a long-term vision of ocean management.


Introduction
Oceanographic connectivity can be broadly defined as the existence of "a specific path between identified places" [1].This definition relies heavily on the hydrodynamics that govern the transport of elements in the oceans.Oceanographic connectivity is associated with inert elements (e.g., pollutants, microplastics), and living elements (e.g., fish larvae), which are exposed to seasonal variability and depth-related variations in the hydrodynamic fields.Many elements are only transported by the marine currents, with some of them having a limited control over their transportation like early life stages of many marine species.Other elements have complete control over the currents and dominate their dispersal in the oceans like, for example, highly migratory species.For marine beings, connectivity involves and impacts diverse biological scales (i.e., individuals, species, populations, communities, and ecosystems).Connectivity is fundamental to guarantee gene flow and the spatial structure and dynamics of populations within the network and subject to extreme oceanic events [2,3].The genetic mixing brought by the newly arrived individuals contributes to the stability and resilience of the population against environmental shifts and pathogen aggressions [4].Besides, marine connectivity events may produce impactful changes in the ecological status of areas, such as massive arrivals of invasive species (e.g., the Pacific oyster in the Wadden Sea [5]).
Estimates of connectivity can be provided from particle transport modelling, a numerical tool in which oceanography is the indispensable input.Particle transport modelling gives substantial information on the complexity of interaction between particles and its environment.This type of model is a powerful, well-used, and adaptable method to estimate connections among marine areas, and thus, representing connectivity [6,7].They are desk-based and allow the simulation of a significant number of transported particles.These tools have often been used for the transportation of small elements with little or no ability to move, which interfere with current advection.For connectivity studies, they have provided relatively accurate estimates of trajectories between locations in the oceans and the flux of individuals within trajectories.In a research context seeking for sustainability of human activities and predicting the effects of global warming, many modelling studies have examined the connectivity of marine populations, through the larval dispersal of commercial species or umbrella species (e.g., habitat formation) [6,7].Results from these models have contributed to understand marine connectivity and provided information for marine spatial conservation planning such as the prioritisation of zones in Marine Protect Area (MPA) management [8].
Across time, marine connectivity is likely to change due to the high temporal variability in the atmospheric conditions forcing surface ocean circulation [9].The main circulation current is bound to have variation in its velocity amplitudes, its geographical ranges (e.g., the Northern Current in the Mediterranean Sea [10]), and its directions over time, from a day-by-day basis to decades (e.g., the Gulf Stream [11]).Mesoscale circulations, such as gyres and fronts, have similar fluctuations in addition to their occurrence in the regional hydrodynamics.Circulation patterns may be accentuated and modified by climate change, to which, observations and models have shown and projected perturbations [12][13][14].Temporal variability in water properties (e.g., temperature) affects ecological factors (e.g., reproductive dynamics and larval pelagic duration) and therefore, affects marine connectivity of populations, as has been shown by many biophysical modelling studies [15][16][17].
Knowledge of the seasonal and interannual variability in connectivity processes sheds light on the sites that have stable or intermittent connections over time [18,19] and is useful for management and monitoring marine environments at middle and long term temporal scales.Connectivity estimates are presented or addressed on a temporal unit basis (e.g., comparison of connectivity between years), but the temporal variability in the connectivity over time is less represented.Only a few studies (e.g., [20][21][22] have conducted temporal variability approaches in the context of conservation managing plans.Estimating the variability in connectivity through time is a way to detect the sensitivity of connectivity to time, especially during strong anomalies in marine currents (i.e., years with strong events of deep-water formation, cascading, heatwaves, and storms) and during upwelling processes.Additionally, it is also a way to identify zones where connectivity can be similar among zones or where the patterns look unique.Overall, information on the temporal variability in connectivity combined with ecological traits from selected marine species can support the policy-makers with conservation decisions.Nonetheless, it lacks of concise and adapted indicators.
In this study, we aimed to define three theoretical indicators for characterising the temporal variability in occurrence, flux, and frequency of connectivity in an area.A first calculation of the indicators was set up in semi-theoretical context using the estimates of oceanographic connectivity in the northwestern (NW) Mediterranean Sea across time obtained by numerical modelling (i.e., coupling a particle transport model with regional hydrodynamics).In addition, we investigated the impact of two relevant features impacting connectivity on the indicators estimates: the transport duration and the depth of particle sources.The proposed three indicators are thought to be ecological metrics for assessing the variability in the connectivity for zones of spatial conservation interest, and for enhancing and improving the monitoring and governance of marine resources.

Definition of connectivity indicators
We defined three indicators to characterise the temporal variability in connectivity at a location: occurrence, flux, and frequency (Table 1).Their calculations relied on the creation of links (i.e., a connection) established by at least one transported element, here a particle, between source and destination for a transport duration D. Application of the indicators focused solely on either the source or destination perspective.Links were stored in a dataset Table 1.Definition of the three connectivity indicators at a focus site i and its linked sites j.

Indicatorª Definition Notation and calculation Range
Occurrence P X>1 : The proportion of repeated links over the cumulated in-degree/out-degree* at focus site n link : The time-average of cumulated in-degree/out-degree* at focus site Flux β H q : Flow rate of particle over time (beta Hill number with q = 5, S1 Fig) P N|X = 1 : Proportion of particles flowing in non-repeated links over the total number of connecting particles at a focus site Frequency F: Over a continuous period of particle release, the frequency of linking M link : Over a continuous period of particle release, the frequency of maximum uninterrupted link duration X ij is the number of occurrence of a repeated connection between focus sites i and linked site j, f i is the total number of linked sites with focus site i, T i is the number of release times with established linked from focus site i, S Xij a binary index indicating the repetition of links between focus and linked sites across time, d is the total number of release times, N i,t is the out-strength/in-strength from the focus site i at release time t, N ij is the link weight between focus site i and linked site j across time, N ij|t is the link weight from a focus site i to a linked site j at release time t, N i is the cumulated out-strength/in-strength across time, and x ij is a vector of values representing the durations of uninterrupted links between a focus site i and a linked site j.
that contained information on the Source ID, the Destination ID, the transport duration D, the particle release depth (binarized as near-bottom or near-surface), the particle release time, and the considered time units (month or year).We did not consider accumulation of links that were established before the time D in the calculations.Calculations of indicators did not treat 'retention' links differently from other links (i.e., when a source was connected to itself).Finally, for the Frequency indicator, a constant time step between releases in the particle transport was needed.
The three indicators ranged within specific values, generally between 0 and 1 (see Table 1), and some values being influenced by others (e.g., P x>1 to n link , and P N|X = 1 to P x>1 , and M link to F). Occurrence indicators highlighted whether the focus sites had numerous connections (n link ) and whether the connections were repeated (P x>1 ) across time.The occurrence connectivity indicators n link and P X>1 relied on the degree measure from the Graph theoretic approach [23], considering the source perspective (out-degree) and the destination perspective (in-degree).If n link = f, the focus site was connected to all linked sites every time.On the other hand, if n link = 0, the focus site never linked to another site.Under the condition that n link > 0, P x>1 = 0 would mean that links were established once across time, while P x>1 = 1 would mean that connections were established every time between the focus and linked sites.Flux indicators represented the variations of particle quantity flowing through the links ( β H 5 ) and the impact of particle quantity in unrepeated links (P N|X = 1 ).We stressed out that calculation of the N variables relied on strength and link weight measures from the graph-theoretic approach [23].The out-strength and in-strength measures are used from a source and destination perspectives, respectively.For the Flux connectivity indicator β H 5 , we adapted the application of Beta Hill number β H q , that usually measured the variation of species diversity and abundance within samples [24].Hill numbers have been found useful for marine management [25].We selected the Hill order q = 5 to emphasise the sensitivity of the Hill number to variations in links with an abundant flux.This value was chosen using an elbow approach during the analysis of β H variations with q (S1 Fig) based on particle transport simulations described in the next section.If β H 5 = 1, the particle quantity in the link scarcely changed across time, but if β H 5 approached T, the variation of particle quantity significantly changed.Focusing on the impact of non-repeated link, we interpreted that if P N|X = 1 = 1, all particles flowed in the nonrepeated links and if P N|X = 1 close to 0, then, a majority of the particles flowed in the repeated links.Frequency indicators highlighted the consistency of link establishment over time, providing insight into the duration of links throughout the studied period.Frequency indicators are calculated exclusively when there is a constant time step in particle releases.If F = 1, connections were repeated across all linked sites over time.When M link = 1, the links remain continuously established throughout the entire time period.Smaller values of M link decrease the likelihood of connections being sustained between release times.

Biophysical modelling
Context of modelling.The indicators were calculated in a semi-theoretical example using the results of a Particle Transport Model (PTM) set up in the NW Mediterranean Sea (Fig 1).Such context was configured to challenge the indicators against complex ocean hydrodynamics.In the NW Mediterranean Sea, there is a well-studied southward current following the Spanish Mediterranean coastlines [26].Mesoscale circulations introduce variability in the local circulation (e.g., gyre position, sizes, meanders along the general current), make the general circulation pattern spatially fluctuate over time [27], and affect connectivity [15].
To apply the indicators that we defined (Table 1), the present case study restricted their calculations on selected areas out of 34 zones (see Fig 1) in the NW Mediterranean Sea.These 34 zones officially corresponded to marine protected areas [28,29], but their management status was discarded from the models, and the mention of these MPAs was kept to "zones".More particularly, four zones were selected for particle releases in the particle transport models (PTMs).These four zones have similar depth distribution (300-600 m) along the continental slope and sided a well-known hydrodynamic front, having an interannual impact on larval connectivity [15].
Hydrodynamics.To force the advection of particles, the PTM used three-dimension velocity fields of a climatological (i.e., average circulation) and an interannual Regional Ocean Modelling Systems (ROMS [30]) run for the NW Mediterranean Sea within the domain 38N -43.69˚N, 0.65˚W-6.08˚(Fig2).The velocity fields from the hydrodynamic model are provided daily for the climatological year and between the years 2006 and 2020 on a grid with a spatial resolution around 2 km and vertically discretized over 40 sigma layers.Details on the implementation and validation of the climatological year are provided in [31,32].The 2006-2020 model has been forced with interannual atmospheric data provided by the Copernicus' Climate reanalysis database (ERA-Interim from 2006 to 2016 and ERA-5 from 2016 to 2020 [33,34]) and provided by the Cross-Calibrated Multi-Platform (version 3) from the Remote Sensing System [35].Data were standardized to the ROMS model grid and to a time dimension of 6 hours for the air temperature, cloud cover, freshwater flux, wind components, and surface pressure and 12 hours for the shortwave and longwave radiations, and precipitation.Further details on the ROMS calibration and validation are provided in [15] For the implementation of the PTM with interannual variability, we used current velocity fields, whose years were selected after an Empirical Orthogonal Function (EOF) analysis on the yearly surface Mean Kinetic Energy (MKE).MKE was transformed, such as MKE values for the EOF analysis were squared-root cosine-weighted standardised anomalies (1): MKEðX;Y;tÞÀ � MKEðtÞ s MKEðtÞ � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi cosðY � p=180Þ p , with � MKEðtÞ and s MKEðtÞ , the spatial average and standard deviation of MKE, respectively, X and Y, the longitude and latitude in degree, respectively, and t, the year between 2006 and 2018.
EOF is a spatiotemporal mathematical method identifying patterns of variability in oceanographic data [36] and relies on the principal component analysis.The two first EOF mode for MKE (explaining 38.5% of the total variance) revealed two main spatial structures (Fig 2).The first mode of the EOF (EOF1) represented a meander structure, which is associated with the dynamics of the Northern current, flowing north-easterly over the Gulf of Lion.The second Transport modelling.Particle transport was modelled with the open source Opendrift software [37], using a Runge-Kutta 4 th order scheme for particle advection.The time step for trajectory calculation was set to one hour and the particle positions were saved daily.If a particle hit land, the particle was deactivated and stayed at the same position until the end of the simulation.
Passive particles were randomly released within the spatial delimitation of the four selected zones 2, 8, 18, and 21 (Fig 1).For simulations carried out in the selected years 2008, 2014, and 2018, the depth of the releases was set near the surface (5 m below the surface), while for simulations carried out with the climatological year 'clim', the depth of the releases was set near the surface and near the bottom (5 m above the bottom).The simulation of the particles' transport started every first day of the month in the three separate years.The representativeness of the particles as early life stages of marine species has been limited when setting up the models to keep the approach as theoretical as possible.However, the transport duration lasted a maximum of 30 days, which was a value encompassing the pelagic duration of many marine species (e.g., [38]).Releases at higher frequency (e.g., every day of a month) were not implemented as we were not interested in the temporal variability in the connectivity at small time scales.A total number of 88,763 particles were released each month.This number represented more

Calculations of connectivity indicators
The connectivity indicators were calculated using different combinations of particle simulations to isolate effects of temporal variability and PTM parametrisations (Table 2).Except for the analysis of the transport duration effect on the indicators (see below), we selected the time D of 10 days for calculating the indicators.In the present study, D represented the moment a particle stopped being transported, which, in larval connectivity studies, corresponds to either the moment larvae would settle after reaching competent stages and potentially favourable habitats or the moment they would die.The settlement age is dependent to the species and to environmental contexts, hence several values of D were used in a sensitive analysis.Our study used three discontinuous years representing extreme conditions (2008, 2014, and 2018), which did not respect an assumption of consistency in the time steps between releases.Accordingly, we calculated Frequency indicators (see Table 2) using monthly releases spanning three years (2008, 2014, and 2018).
In addition to focusing on the connection from a source perspective, we also calculated the indicators described in Table 1 from the reception perspective to illustrate their bidirectional use (Table 2).The calculation remained the same, but focus sites were destination zones and linked sites were source zone under the reception perspective.Connectivity indicators under the reception perspective sum up which destinations are the most and least connected to sources.
We tested the PTM parametrisation by calculating and comparing indicators for a selected transport duration D (5, 10, 15, and 30 days) and a binary set of release depths (near the sea bottom and near the surface).Only two release depths were investigated following the rationales of making release in two opposed sets of hydrodynamic fields with currents generally being slower and steadier close to the bottom than near the surface.

Results
The connectivity between the four sources and 35 destinations are summarised in connectivity matrices (Fig 4 ) for the three years, an average of the three years ('Avg'), and the climatological year ('clim').The climatological year 'clim' approached connectivity patterns that were provided by particle transport simulation in an average circulation field of the NW Mediterranean Sea.For D = 10 days, we found that: first, for the three years, the sources did not connect with all the 34 destinations (8 zones with the IDs 25-27 and 30-34 were not connected to the sources).Second, the connectivity from the sources could be limited to a range of destinations.For example, source zones 2 and 8 were not connected to southern destinations (zones identified with an ID � 17).Lastly, on a yearly basis (i.e., 2008, 2014, and 2018), the connectivity matrix also highlighted the connectivity variability per sources.For example, source zone 2 was well connected throughout the three selected years while the connectivity of source zone 21 is low and limited to up to 8 zones in year 2014.This diversity of temporal variability in connectivity is quite distinct from the average of connectivity across time or from connectivity estimated with the hydrodynamic in a climatological.In contrast, when examining a 3-year average of connectivity, it suggested that source zone 18 had connections with 21 zones.However, upon analysing the results derived from the climatological year, we only counted two connections between zone 18 and destination zones.

Indicators and temporal variability
To provide a clearer overview of the indicators and differences among the focus sites, we plotted indicators of the same category on a plane.Source perspective.The four source zones of particles presented diverse temporal variability in connectivity, as shown by the scattered distribution of the source's indicator values (Figs 5-7).
The Occurrence indicators (Fig 5 ) showed that connectivity among zones in the NW Mediterranean Sea is highly diverse.For instance, source zone 21 tended to connect a single time to a destination more often than repeatedly (P X>1 = 0.46) over the three simulated years and, on average, was linked to fewer destinations (n link = 4.3) than the other sources.Conversely, source zone 2 was linked at least twice to the same destinations (P X>1 = 1) and connected to more destinations (n link = 14) than the other source zones.The same can be observed by analysing the connectivity matrix (see Fig 4).Source zone 21 connected with 7 destinations once, mostly in the year 2014 (5 single links) and source zone 2 connected at least twice with 16 destinations (thrice with 10 destinations).These Occurrence indicators alone already highlight the relevance of zone 2 position for likely creating persistent links with others.On the spatial dimension, zone 2 is one of the northernmost sites of the studied area and is geographically close to 11 zones in a radius of 50 km (see Fig 1).Its proximity to several zones can help to have a high n link .However, the proximity of other zones does not necessarily imply n link would always be high.This statement is illustrated by the less connected source zone 8 (n link = 7.67), which was geographically close to 14 zones within a radius of 50 km, while it was located 30 km south of source zone 2.
The Flux indicators (Fig 6 ) revealed the subtleties across time occurring in the established links, and this statement also applies for the links occurring just once.Source zone 21 can be considered a source with limited flow of particles among the four source zones because connections were rare and lowly impactful over the selected years.Indeed, the flow of particles through the single links (P N | X = 1 ) from source zone 21 represented 9.7% of the particle flow.Besides, its flow of particles towards the other destinations was relatively constant ( β H 5 close to 1).These statements were also visible within the connectivity matrix (Fig 4).A total of 125 particles were transported through the single links from source zone 21.To a lesser extent, source zone 18 also had low-impactful single links (P N | X = 1 = 22.1%) but had more variations of particle flows through the repeated links ( β H 5 close to 3).For source zone 2, the Beta Hill number highlighted that the flow of particle relatively and mildly fluctuates between the three selected years.This would mean that the position of zone 2 is not only strategic to connect with several destinations at several times (i.e., high n link ), but also to sustain a consistent connectivity flow across time, making that zone a potential super spreader site.
The Frequency indicators could discern the monthly variability in connectivity in the years 2008, 2014, and 2018 and among the four sources (Fig 7).On the first hand, the monthly connectivity in the years 2008 and 2014 had opposite patterns.The low values of the indicators in the year 2008 indicated a less frequent establishment of links along the 12-month perisod (F 0.1) and a lower probability of consecutive connections (M Link = 0.094) than for the other years.It was the opposite for the year 2014: links lasted a maximum of 1.5 months (M Link = 0.13) and occurred 2.1 times in 12 months (F ~0.18).Taken separately (S3 Fig), the matrices of connectivity also represent these distinct frequency features of connectivity.On the other hand, the Frequency indicators revealed that source zones have relatively comparable variation of connectivity through 12 months.Indeed, in the three studied years, the Frequency indicators of source zones 18 and 21 were generally lower (F < 0.14 and M Link < 0.11) in comparison to the source zones 2 and 8 (F > 0.14 and M Link > 0.12).
Reception perspective.An alternative perspective was to consider the destination of particles at the zones of interest (i.e., 2, 8, 18, and 21) as focus sites (Fig 8).We found that all destination zones, except zone 21, had repeated connections (P X = 1 > 0.5) and are bound to receive a similar flow of particles every time they are connected ( β H 5 2 [1-1.15]).Particular cases of temporal variability in connectivity were highlighted such as the high proportion of

Sensitivity of the indicators to model parametrisation
Release depths.The depth of releases impacted the indicators in three different ways (Fig 9).First, null or minor changes in the connectivity indicators across time and sources (e.g., n link ) were detected in the indicators.The Occurrence indicator n link shows that the number of links is likely to stay the same if the dispersal occurred near bottom or near surface (e.g., for source zone 18, n link = 1.29 with surface simulations and n link = 1.75 with bottom simulations).While the changes in n link remained small with the two depth simulations, the destinations of particles can be subject to modification.For example, source zone 18 connected to zones 21 and 22 with bottom simulations and to zones 16 and 18 with surface simulation (S4 Fig).Second, the indicators for a few source zones were modified (e.g., β H 5 and P N|X = 1 for zones 2 and 21 in Fig 9B ).For example, the releases made near-bottom decreased and increased β H 5 in source zones 2 and 21, respectively, in comparison to releases made near-surface.Finally, the depth of release marked a clear dichotomy in indicators for most (if not all) zones (e.g., P X > 1 and M link in Fig 9A and 9C).The probability that the source had a repeated link to a destination P X>1 decreased by 0.17 for releases near the bottom (P X>1 = 1.00) in comparison to releases near the surface (P X>1 = 0.83 on average) for all sources except source zone 21.The dichotomy was also well seen in the Frequency indicator M Link , with longer continuous connections for particle released near the bottom (M Link = 0.55 on average) than near the surface (M Link = 0.19 on average).
Transport duration.Analysing the sensitivity of the connectivity indicators to the transport duration through the selection of times D revealed that there are multiple possibilities of connectivity variations depending on the source ( Fig 10).First, the indicator value increased/ decreased with increasing D from 5 to 30 days.For example, P X>1 from 0.93 to 0.66 in source zone 2 and n link from 3. 0 to 6.0 in source zone 8. Second, some indicators had small or big variations in their values without linear relationship to the transport duration, e.g., source zone 8 with β H 5 2 [1.0, 2.18] and source zone 21 with P N|X = 1 2 [0.01, 0.54].Third, the zones could also be grouped according to relatively low and high amplitudes of indicators, suggesting that there was a different impact of the transport duration based on the source zone and their geographical position in the NW Mediterranean Sea.For example, the n link was lower for southern source zones 18 and 21 (i.e., on average 1.5 and 1.6, respectively) than for northern source zones 2 and 8 (i.e., on average 3.6 and 4.3, respectively).

Discussion
The proposed indicators allowed us to characterise the variability in connectivity from sources and destinations and contributed to identifying zones for spatial conservation management that have potentially high oceanographic connectivity through time.

Temporal variability in the connectivity
The temporal variability in the connectivity was influenced by the spatial distribution of zones along the NW Mediterranean margin.On the one hand, their geographical distribution in the area played a role in the variability of connectivity due to their exposure to the dominant water circulation, as often observed in other regions [40,41].The formation of multiple connections over time, defined by the Occurrence indicator P X>1 and its inverse (i.e., 1-P X>1, when a source zone tends to connect a single time to another zone) highlighted how some zones were more exposed to the temporal variability in currents than other zones.For example, single connections may occur when particles are advected in velocity fields (by their directions or amplitudes) established under abnormal circumstances (e.g., under the influence of a cyclone [42]).A low P X>1 implied that the water circulation would normally isolate one zone from the rest and that, for some years, anomalies in the circulation pattern would suspend this situation and allow connectivity (e.g., results on zone 21).In our case, the anomaly was due to a strong southward current in 2014 contrary to the circular current observed in the years 2008 and 2018.Such results recall the importance of having hydrodynamic models that well represent the mesoscale structures of the circulation.A fine resolution of the model is a prerequisite for a good representation of connectivity [43], especially when larval dispersal occurred in coastal waters.Our study was less sensitive to the hydrodynamic model resolution of 2 km because it was conducted in open waters and for zones above the continental slope.However, if a finer model were available, the differences in connectivity indicators between zones would probably be more pronounced.On the other hand, the distance between the zones (i.e., the spacing between sources and destinations) has a mild impact on the connectivity temporal variability: a network of close zones did not necessarily imply a higher occurrence of connections (e.g., low n link for source zone 8) in our theoretical approach.Spacing is a relevant criterion for defining zones as MPAs, following proposed recommendations [44,45].Yet, if the geographical distribution of spatially close zones always overlaps with interannually variable mesoscale structures, as in zone 8, connectivity and its variability may be restricted and less relevant for MPA performance.From this theoretical finding, we suggest that the persistence of local hydrodynamics should also be investigated when designing spatially-close MPAs.
The present indicators can be used to characterise zones on the basis of their success in establishing connections across time.For instance, a source zone with numerous connections (high n link ) that are repeated over a unit of time (high F) can be defined as a 'persistent' super spreader, using terminology of the graph-theoretic approach.If it has only one repeated connection (low n link ) to a destination, it can be defined as an 'persistent' restricted spreader.In terms of the reception perspective, a destination zone that maintained several connections (high n link ) from time to time (low F) can be called an 'occasional' super-receiver.Moreover, by examining indicators from both source and destination perspectives, we can identify the presence of "persistent" loops of connectivity across time.A recurring loop pattern would be favourable to the population sustainability and align with predictions from metapopulation models [46].The classification of zones based on characteristics of temporal variability in connectivity has a practical outcome for stakeholders' understandings, especially when designing a network of MPAs whose implementation is long term.The current study presented these indicators in specific and limited conditions, consequently we are not able yet to provide thresholds in demand by policy-makers [47].The indicators should be stressed and tested by multiple applications in the future in order to define threshold values that are either robust and relevant wherever and whenever they are used or specific to the case study (e.g., F MSY [48]).As for now, these indicators enhanced the need for further analysis of certain zones to understand whether they are relevant and essential for regional connectivity.In fact, the indicators could reveal unusual connectivity for certain zones (e.g., zone 18) where it might be strategic or not to make conservation decisions.For instance, would a protected area that tends to be connected only once and by a few particles to a non-isolated MPA destination (i.e., one that is already well connected to other MPA sources) be a priority for implementation or modification?The answer may be complex, as the impact of the presence of an MPA on connectivity in a network may be quite significant [49].Calculating the temporal variability in connectivity can provide another layer of criteria to inform management decisions when previous assessments are inconclusive.

Connectivity indicators with ecological implications
The calculations of the indicators were based on a semi-theoretical approach and have yet to be applied to larval transport studies with specific species characterisations.The PTM parametrisations approximated key ecological traits of early-life stages of marine species, affecting connectivity: 1) depth of particle release as the spawning depth of demersal and pelagic species and also surfacing vertical larval behaviour, and 2) transport duration as the pelagic life duration (PLD) of early life stages (including the different forms of the pelagic stages, such as eggs or larvae).The sensitivity of the indicators to the parameters showed interesting results for using indicators to compare the connectivity of marine species and to understand the impact of global warming on connectivity.Firstly, the indicators are likely to be helpful in showing when early life stages with short or long PLD, or with a benthic or epipelagic habitat would be affected in an environment with temporal variability.The expected results would be in line with the many existing studies addressing the issue of ecology in connectivity [38,50,51].In addition, these indicators shall assess whether a zone allows connectivity for multiple species or for a single species once ecological traits are included in the transport models.It implies divergent or convergent connectivity patterns as function of the species life habits [51,52].This would permit to understand if the MPA network is ecologically coherent [53], meaning that it is adequately designed for various groups of species.Secondly, by measuring and comparing connectivity indicators over two periods separated by a known regime shift, the temporal variability in the connectivity due to climate change can be estimated.In the NW Mediterranean Sea, an increase in the water temperature has been observed over the last two decades [54] and correlated with various changes in marine community in the Mediterranean region [55][56][57].In the transport of early life stages, the time spent in the pelagic waters could shorten or lengthen according to species-dependent metabolic and behavioural traits associated with water temperature [58].Consequently, connectivity is likely to change over time.
Further work on the calculation of indicators should help to improve the measurement of the temporal variability in connectivity when including the ecological dimension.The context of the study limited this investigation but the development and calculation of indicators would be welcome.Firstly, we lacked information on the temporal quality of connectivity, which can nevertheless be approached within the calculation of the Flux indicator.The quality of the connectivity should take into account the loss of particles before, during and after the transport, and is strongly, if not exclusively, related to the intrinsic ecology of the species (e.g., spawning success, mortality, growth, or reproduction rates).These ecological traits are essential information for determining with greater precision the temporal variability in connectivity of a given set of populations in a geographic network of MPAs [59].Second, a useful Occurrence indicator would be to indicate whether an MPA is receiving and sourcing dispersing particles.This metric would make it possible to understand whether an MPA is receiving more/equal/less larvae than it is delivering to other MPAs or whether it is acting as a recurring stepping stone unit over the course of several years [60].In fact, the methodology of the graph-theoretic approach [61], already used to design and evaluate MPA networks [48,60], should help to develop indicator calculations that focus on the temporal variability in connectivity.Third, the calculation of another Frequency indicator should help to determine the patterns of seasonality in the course of several years and provide meaningful information for the management of temporary MPA closures.At the species level, connectivity can be episodic and seasonal, since marine species have different spawning strategies and life cycles, resulting in peaks of connectivity at certain times of the year.In particular, many species release individuals (i.e., eggs, larvae, spores, propagules) either at a specific time of the year (e.g., winter in the case of Nephrops norvegicus in the NW Mediterranean Sea [62]), relatively frequently throughout the year [63] like the European hake Merluccius merluccius in the NW Mediterranean Sea [64].

Indicators for conservation planning and assessment
The purpose of defining indicators of temporal variability in connectivity is, in part, to support spatial conservation planning.Studies have shown that an optimal MPA network can be estimated from an analysis over a relatively short period of time (i.e., 10 years [65,66]) and that MPA effects are expected to be seen after a few years [67].Such findings suggest that: i) interannual monitoring and assessment of MPAs should be considered, ii) this optimal period of time should be searched for the MPA network in the NW Mediterranean Sea, and iii) the usefulness of the presented connectivity indicators should be questioned by calculating them over the selected optimal time period.
In our study, the zones were 35 deep-sea MPAs in the NW Mediterranean Sea, established for the recovery and conservation of target species of commercial fisheries (i.e., the European hake Merluccius merluccius, the Norway lobster Nephrops norvegicus, the blue and red shrimp Aristeus antennatus) and deep-sea bioengineers (e.g., soft-bodied cold-water corals and gorgonians).Connectivity has not yet been assessed in these MPAs, but monitoring data have begun to be collected and processed [68,69].Calculation of the indicators of temporal variability in connectivity in this network, using a setup of ecological traits in the PTM, would provide benchmarks for MPA assessment and understandings of long-term connectivity in the ocean.They appear to be a good and robust way to compare the connectivity variability of MPAs, notably with the plane representation that visually discriminated MPAs.It is essential to acknowledge that these indicators do not pretend to replace other existing methods of connectivity analysis, e.g., connectivity matrices, graph-theoretical indicator, connectivity portfolio effect [70].In fact, we still need matrices to understand and quantify what is observed with the indicators.And, contrary to graph-theoretic indicators (e.g., betweenness centrality), our indicators do not account for the indirect connections, which occurred when a focus site linked to sites that, in turn, are also connected to others.
Far beyond their use for comparison, indicators can be helpful to redesign MPA networks, adapting their structure to locations with interesting long-term connectivity patterns, and thus, improving the conservation of marine ecosystems.Redesigning MPAs with PTM requires a computational effort, but it helps to optimise the connectivity in a geographic network [71].As long as spatial positions are stored over time, connectivity indicators can be calculated on datasets obtained from other methods involved in the movement tracking, such as Lagrangian drifters (e.g.Argos buoys), which would enhance the aspect of oceanographic connectivity; as well as mark-recaptured tagging of animals, which could help understand the dispersal of the species to different zones for their ecological needs (e.g., feeding, foraging) over years.The application of the indicators is also open to different types of marine spatial management planning, such as the spread of marine invasive species and recruitment to fishing grounds.

Conclusions
Three indicators characterising temporal variability in connectivity have been defined and applied in 35 zones distributed in the NW Mediterranean Sea.These indicators are based on the temporal variability in connectivity established among zones of interest (i.e., MPAs) through particle tracking modelling.They encompass the notion of connectivity Occurrence, Flux and Frequency among zones.These indicators would be relevant for decision-makers during the design process of MPAs or when reporting on their efficiency.Applying these indicators should provide and underline the similarities and differences in terms of connectivity among zones, hence refining a location appropriateness among other candidate zones to management planning.Additionally, a likely use of the indicator comparison is to estimate the capacity of the zone to allow connectivity for species with different ecological traits and the connectivity variation under climate change scenarios (e.g., increasing water temperature).Overall, the indicators are a summary of connectivity variability when multiple time steps are considered in a study, nonetheless thresholds should be defined to better advised the management authorities.We consider that these indicators are applicable irrespectively to the studied regions (here the NW Mediterranean Sea) and to the connected zones (here MPAs), and more importantly, relevant for assessing the MPA network coherency and prioritising management efforts on some MPAs of the network.

Fig 1 .
Fig 1. Zones of particle release and potential destinations in the NW Mediterranean Sea.The network of MPAs is represented by grey polygons.The four zones with labels are the release areas of particles.https://doi.org/10.1371/journal.pone.0297730.g001

Fig 3 .
Fig 3. Principal component (PC) time series of standardized weighted MKE anomalies from the EOF analysis.Values of first PC (A) and second PC (B) are shown for the period 2006-2020.https://doi.org/10.1371/journal.pone.0297730.g003 In 2014, source zone 18 established the majority of connections, totalling 18, in the year 2014, and drove the patterns observed in the 3-year average.These simple examples succinctly demonstrated the purpose of analysing the fluctuations over time in assessing connectivity patterns.

Fig 4 .
Fig 4. Connectivity matrices between the zone sources and destinations across three selected years, their average ('Avg') and a climatological year ('clim').Duration of particle transport D: 10 days.Colour gradient represents the particle abundance flowing from a source to a destination.The grey-shaded boxes stand for zero connected particle.Zone destinations are ordered by their latitude over the NW Mediterranean Sea: Zone 1 being the northernmost and Zone 34 the southernmost.https://doi.org/10.1371/journal.pone.0297730.g004

Fig 5 .
Fig 5. Occurrence indicators n link and P X>1 according to zones as a source of particles.Indicators (see Table 1) represented temporal variability in connectivity across 2008, 2014, and 2018.Duration of particle transport D: 10 days.Focus sites-Zone 2: White, Zone 8: Green, Zone 18: Light brown; Zone 21: Brown.

Fig 8 .
Fig 8. Occurrence (A), Flux (B), and Frequency (C) indicators according to zones as a destination of particles.Indicators of Occurrence and Flux (see Table 1) represented temporal variability in connectivity across 2008, 2014, and 2018.Indicators of Frequency represented temporal variability in connectivity across the 12 months of the year 2014.Duration of particle transport D: 10 days.Focus sites-Zone 2 white, Zone 8: Green, Zone 18: Light brown; Zone 21: Brown.Other zones: Grey.https://doi.org/10.1371/journal.pone.0297730.g008

Fig 9 .
Fig 9. Occurrence (A), Flux (B), and Frequency (C) indicators according to source zones and release depths.Indicators (see Table 1) represented temporal variability in connectivity across 12 months of the climatological year.Duration of particle transport D: 10 days.Focus sites-Zone 2: White, Zone 8: Green, Zone 18: Light brown; Zone 21: Brown.Release depth-near the bottom (round shape) and near the surface (triangle shape).